
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The impact of divorce laws on the equilibrium in the marriage market.
% Reynoso
% April 2024
%
% This file computes the variance-covariance matrix of the estimates and
% the sensitivity moments.
%  
% Data: PSID 1968-1992
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%------------------------ Preliminaries ----------------------------------%

clc
clear all

%--- Indicate location of Replication folder:
%replication_location = 'C:\update_with_your_path';

%--- Working directory
cd ([replication_location,'\Internal_estimation']);

%--- Path for estimates
estimates_dir = [replication_location,'\Outputs\MCD\'];

%--- Path for storing
output_dir = [replication_location,'\Outputs'];

%--- Path for Model
addpath(genpath([replication_location,'\Internal_estimation\Inputs']));

%--- Path for Data
data_dir = [replication_location,'\Data\Inputs'];
mom_imp = importdata([data_dir,'/moments_boot.csv']);
regmom_imp = importdata([data_dir,'/reg_boot.csv']);
moments_boot = [mom_imp.data regmom_imp.data]; 
moments_boot(any(isnan(moments_boot), 2), :) = [];

C=cov(moments_boot); % covariance matrix of bootstrapped data

%--- Load estimation results
%load('outputs_mcd'); 
load([estimates_dir, '\outputs_mcd']);
params=X;

%------------------------ Variace-covariance matrix ----------------------%
%--- definitions
nest=size(params,2);
npar=size(params,2)-n_g*4;
nmom=size(data_moments,1);
[F, model_moments,moment_cond]=moments(n_g, T, lt, dt, H_max, disc,  stoch_earn_f, stoch_earn_m, m_eta_f, m_eta_m, n_eta_f, n_eta_m, home_prod, match_qual, m_th, n_th, gamma, L_max, params, data_moments, data_f, data_m, W);
g=repmat(moment_cond,1,npar); 
increment=0.001.*X; 
h=diag(increment);
X_ptbd_up=repmat(X,nest,1)+h; 
X_ptbd_down=repmat(X,nest,1)-h; 
lambda_rep = repmat(X(npar+1:end),npar,1);
params_ptbd_up=[X_ptbd_up(1:npar,1:npar) lambda_rep];
params_ptbd_down=[X_ptbd_down(1:npar,1:npar) lambda_rep];

%--- Numerical derivative
%// Forward perturbation
g_ptbd_up=zeros(nmom,npar); 
for pt=1:npar
params=params_ptbd_up(pt,:);
[F, model_moments,moment_cond]=moments(n_g, T, lt, dt, H_max, disc,  stoch_earn_f, stoch_earn_m, m_eta_f, m_eta_m, n_eta_f, n_eta_m, home_prod, match_qual, m_th, n_th, gamma, L_max, params, data_moments, data_f, data_m, W);
g_ptbd_up(:,pt)=moment_cond;
end
%// Backward perturbation
g_ptbd_down=zeros(nmom,npar);
for pt=1:npar
params=params_ptbd_down(pt,:);
[F, model_moments,moment_cond]=moments(n_g, T, lt, dt, H_max, disc,  stoch_earn_f, stoch_earn_m, m_eta_f, m_eta_m, n_eta_f, n_eta_m, home_prod, match_qual, m_th, n_th, gamma, L_max, params, data_moments, data_f, data_m, W);
g_ptbd_down(:,pt)=moment_cond;
end
eps_par=increment(1:npar); 
increment_rep=repmat(eps_par,nmom,1); 
deriv_fwd=(g_ptbd_up -g )./(increment_rep);
deriv_bkwd=(g-g_ptbd_down )./(increment_rep);
deriv = (deriv_fwd + deriv_bkwd)*0.5;

%--- Variace-covariance matrix
V=inv(deriv'*W*deriv)*deriv'*W*C*W'*deriv*inv(deriv'*W*deriv);                                % 
diagV=diag(V);
se=sqrt(diag(V));
sensitiv=abs(inv(deriv'*W*deriv)*deriv'*W);

%--- Save
save(fullfile(output_dir,'\MCD\Variance_mcd'),'V', 'deriv', 'sensitiv', 'C', 'se');

%------------------------ Sensitivity moments ----------------------------%

clearvars -except C se sensitiv output_dir
npar=size(se,1);
sd_diag = sqrt(diag(C))';
sd_diag = repmat(sd_diag,npar,1);
sensitiv_norm = sensitiv.*sd_diag;
total2 = sum(sensitiv_norm,2);
rel_sens2 = (sensitiv_norm./total2);
[max2, mom2]=maxk(rel_sens2,5,2); 

sensitivity_moments = zeros(npar,3);

for i=1:npar
   
    %Most important moment
    %M1
    
    if mom2(i,1)>=1 && mom2(i,1)<=24
    sensitivity_moments(i,1)=1;
    end
    
    if mom2(i,1)>=43 && mom2(i,1)<=66
    sensitivity_moments(i,1)=1;
    end
    
    if mom2(i,1)>=85 && mom2(i,1)<=108
    sensitivity_moments(i,1)=1;
    end
    
    if mom2(i,1)>=127 && mom2(i,1)<=150
    sensitivity_moments(i,1)=1;
    end
    
    %M2
        
    if mom2(i,1)>=25 && mom2(i,1)<=33
    sensitivity_moments(i,1)=2;
    end
    
    if mom2(i,1)>=67 && mom2(i,1)<=75
    sensitivity_moments(i,1)=2;
    end
    
    if mom2(i,1)>=109 && mom2(i,1)<=117
    sensitivity_moments(i,1)=2;
    end
    
    if mom2(i,1)>=151 && mom2(i,1)<=159
    sensitivity_moments(i,1)=2;
    end
    
    %M3
        
    if mom2(i,1)>=34 && mom2(i,1)<=42
    sensitivity_moments(i,1)=3;
    end
    
    if mom2(i,1)>=76 && mom2(i,1)<=84
    sensitivity_moments(i,1)=3;
    end
    
    if mom2(i,1)>=118 && mom2(i,1)<=126
    sensitivity_moments(i,1)=3;
    end
    
    if mom2(i,1)>=160 && mom2(i,1)<=168
    sensitivity_moments(i,1)=3;
    end
    
    %M4
    
    if mom2(i,1)>=169 && mom2(i,1)<=252
    sensitivity_moments(i,1)=4;
    end
    
    %Second most important moment
    %M1
    
    if mom2(i,2)>=1 && mom2(i,2)<=24
    sensitivity_moments(i,2)=1;
    end
    
    if mom2(i,2)>=43 && mom2(i,2)<=66
    sensitivity_moments(i,2)=1;
    end
    
    if mom2(i,2)>=85 && mom2(i,2)<=108
    sensitivity_moments(i,2)=1;
    end
    
    if mom2(i,2)>=127 && mom2(i,2)<=150
    sensitivity_moments(i,2)=1;
    end
    
    %M2
        
    if mom2(i,2)>=25 && mom2(i,2)<=33
    sensitivity_moments(i,2)=2;
    end
    
    if mom2(i,2)>=67 && mom2(i,2)<=75
    sensitivity_moments(i,2)=2;
    end
    
    if mom2(i,2)>=109 && mom2(i,2)<=117
    sensitivity_moments(i,2)=2;
    end
    
    if mom2(i,2)>=151 && mom2(i,2)<=159
    sensitivity_moments(i,2)=2;
    end
    
    %M3
        
    if mom2(i,2)>=34 && mom2(i,2)<=42
    sensitivity_moments(i,2)=3;
    end
    
    if mom2(i,2)>=76 && mom2(i,2)<=84
    sensitivity_moments(i,2)=3;
    end
    
    if mom2(i,2)>=118 && mom2(i,2)<=126
    sensitivity_moments(i,2)=3;
    end
    
    if mom2(i,2)>=160 && mom2(i,2)<=168
    sensitivity_moments(i,2)=3;
    end
    
    %M4
    
    if mom2(i,2)>=169 && mom2(i,2)<=252
    sensitivity_moments(i,2)=4;
    end
    
    %Third most important moment 
    %M1
    
   if mom2(i,3)>=1 && mom2(i,3)<=24
    sensitivity_moments(i,3)=1;
    end
    
    if mom2(i,3)>=43 && mom2(i,3)<=66
    sensitivity_moments(i,3)=1;
    end
    
    if mom2(i,3)>=85 && mom2(i,3)<=108
    sensitivity_moments(i,3)=1;
    end
    
    if mom2(i,3)>=127 && mom2(i,3)<=150
    sensitivity_moments(i,3)=1;
    end
    
    %M2
        
    if mom2(i,3)>=25 && mom2(i,3)<=33
    sensitivity_moments(i,3)=2;
    end
    
    if mom2(i,3)>=67 && mom2(i,3)<=75
    sensitivity_moments(i,3)=2;
    end
    
    if mom2(i,3)>=109 && mom2(i,3)<=117
    sensitivity_moments(i,3)=2;
    end
    
    if mom2(i,3)>=151 && mom2(i,3)<=159
    sensitivity_moments(i,3)=2;
    end
    
    %M3
        
    if mom2(i,3)>=34 && mom2(i,3)<=42
    sensitivity_moments(i,3)=3;
    end
    
    if mom2(i,3)>=76 && mom2(i,3)<=84
    sensitivity_moments(i,3)=3;
    end
    
    if mom2(i,3)>=118 && mom2(i,3)<=126
    sensitivity_moments(i,3)=3;
    end
    
    if mom2(i,3)>=160 && mom2(i,3)<=168
    sensitivity_moments(i,3)=3;
    end
    
    %M4
    
    if mom2(i,3)>=169 && mom2(i,3)<=252
    sensitivity_moments(i,3)=4;
    end
    
end

sensitivity_moments = cell(npar,3);

for i=1:npar
    
    %Most important moment
    %M1
    
    if mom2(i,1)>=1 && mom2(i,1)<=24
    sensitivity_moments(i,1)={'mp'};
    end
    
    if mom2(i,1)>=43 && mom2(i,1)<=66
    sensitivity_moments(i,1)={'mp'};
    end
    
    if mom2(i,1)>=85 && mom2(i,1)<=108
    sensitivity_moments(i,1)={'mp'};
    end
    
    if mom2(i,1)>=127 && mom2(i,1)<=150
    sensitivity_moments(i,1)={'mp'};
    end
    
    %M2
        
    if mom2(i,1)>=25 && mom2(i,1)<=33
    sensitivity_moments(i,1)={'sah'};
    end
    
    if mom2(i,1)>=67 && mom2(i,1)<=75
    sensitivity_moments(i,1)={'sah'};
    end
    
    if mom2(i,1)>=109 && mom2(i,1)<=117
    sensitivity_moments(i,1)={'sah'};
    end
    
    if mom2(i,1)>=151 && mom2(i,1)<=159
    sensitivity_moments(i,1)={'sah'};
    end
    
    %M3
        
    if mom2(i,1)>=34 && mom2(i,1)<=42
    sensitivity_moments(i,1)={'dp'};
    end
    
    if mom2(i,1)>=76 && mom2(i,1)<=84
    sensitivity_moments(i,1)={'dp'};
    end
    
    if mom2(i,1)>=118 && mom2(i,1)<=126
    sensitivity_moments(i,1)={'dp'};
    end
    
    if mom2(i,1)>=160 && mom2(i,1)<=168
    sensitivity_moments(i,1)={'dp'};
    end
    
    %M4
    
    if mom2(i,1)>=169 && mom2(i,1)<=252
    sensitivity_moments(i,1)={'earn'};
    end
    
    %Second most important moment
    %M1
    
    if mom2(i,2)>=1 && mom2(i,2)<=24
    sensitivity_moments(i,2)={'mp'};
    end
    
    if mom2(i,2)>=43 && mom2(i,2)<=66
    sensitivity_moments(i,2)={'mp'};
    end
    
    if mom2(i,2)>=85 && mom2(i,2)<=108
    sensitivity_moments(i,2)={'mp'};
    end
    
    if mom2(i,2)>=127 && mom2(i,2)<=150
    sensitivity_moments(i,2)={'mp'};
    end
    
    %M2
        
    if mom2(i,2)>=25 && mom2(i,2)<=33
    sensitivity_moments(i,2)={'sah'};
    end
    
    if mom2(i,2)>=67 && mom2(i,2)<=75
    sensitivity_moments(i,2)={'sah'};
    end
    
    if mom2(i,2)>=109 && mom2(i,2)<=117
    sensitivity_moments(i,2)={'sah'};
    end
    
    if mom2(i,2)>=151 && mom2(i,2)<=159
    sensitivity_moments(i,2)={'sah'};
    end
    
    %M3
        
    if mom2(i,2)>=34 && mom2(i,2)<=42
    sensitivity_moments(i,2)={'dp'};
    end
    
    if mom2(i,2)>=76 && mom2(i,2)<=84
    sensitivity_moments(i,2)={'dp'};
    end
    
    if mom2(i,2)>=118 && mom2(i,2)<=126
    sensitivity_moments(i,2)={'dp'};
    end
    
    if mom2(i,2)>=160 && mom2(i,2)<=168
    sensitivity_moments(i,2)={'dp'};
    end
    
    %M4
    
    if mom2(i,2)>=169 && mom2(i,2)<=252
    sensitivity_moments(i,2)={'earn'};
    end
    
    
    %Third most important moment 
    %M1
    
    if mom2(i,3)>=1 && mom2(i,3)<=24
    sensitivity_moments(i,3)={'mp'};
    end
    
    if mom2(i,3)>=43 && mom2(i,3)<=66
    sensitivity_moments(i,3)={'mp'};
    end
    
    if mom2(i,3)>=85 && mom2(i,3)<=108
    sensitivity_moments(i,3)={'mp'};
    end
    
    if mom2(i,3)>=127 && mom2(i,3)<=150
    sensitivity_moments(i,3)={'mp'};
    end
    
    %M2
        
    if mom2(i,3)>=25 && mom2(i,3)<=33
    sensitivity_moments(i,3)={'sah'};
    end
    
    if mom2(i,3)>=67 && mom2(i,3)<=75
    sensitivity_moments(i,3)={'sah'};
    end
    
    if mom2(i,3)>=109 && mom2(i,3)<=117
    sensitivity_moments(i,3)={'sah'};
    end
    
    if mom2(i,3)>=151 && mom2(i,3)<=159
    sensitivity_moments(i,3)={'sah'};
    end
    
    %M3
        
    if mom2(i,3)>=34 && mom2(i,3)<=42
    sensitivity_moments(i,3)={'dp'};
    end
    
    if mom2(i,3)>=76 && mom2(i,3)<=84
    sensitivity_moments(i,3)={'dp'};
    end
    
    if mom2(i,3)>=118 && mom2(i,3)<=126
    sensitivity_moments(i,3)={'dp'};
    end
    
    if mom2(i,3)>=160 && mom2(i,3)<=168
    sensitivity_moments(i,3)={'dp'};
    end
    
    %M4
    
    if mom2(i,3)>=169 && mom2(i,3)<=252
    sensitivity_moments(i,3)={'earn'};
    end
    
end

save(fullfile(output_dir,'\MCD\Sensitivity_mcd'),'sensitivity_moments');









